Tarea 1
Tarea 1: Foundations
Integrantes :
- Nicolás García Ríos
Cuerpo Docente:
- Profesor: Felipe Bravo M.
- Auxiliar: Sebastian Bustos e Ignacio Meza D.
Fecha límite de entrega:
Índice:
Objetivo
Bienvenid@s a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en análisis exploratorio de datos y conceptos introductorios de probabilidades. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fín de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Instrucciones:
- La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
- La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
- El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
- Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
- No serán revisadas tareas desarrolladas en Python.
- Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
- Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.
Primera Parte: Preguntas Teóricas
A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas preguntas de forma breve, no más de 4 o 5 lineas.
Pregunta 1:
¿Por qué la estadística es importante?, ¿Que nos permite realizar con los datos?. De algún ejemplo.
Porque la estadística es una herramienta que trabaja el análisis e interpretación de datos, permitiendonos Describir, Decidir y Predecir el mundo que nos rodea.
Un ejemplo visto en clase es el estudio de las grasas saturadas. Una idea común es que comer muchas grasas saturadas, o incluso consumir un poco de ellas es muy perjudicial a la salud. Pero resulta que después de un estudio estadístico se pudo determinar que no solo no son malas, si no que hasta cierto punto estas son beneficiosas para la salud.
De no ser por la estadística, no podríamos haber descrito realmente la influencia de las grasas saturadas sobre la salud, y nos hubieramos basado en un pensamiento más Eurístico que Empírico para la desicion de dietas, nutrición, etc.
Pregunta 2:
Un amigue cercano a usted le comenta que le preocupa salir a la calle cuando hay ofertas en los helados, esto debido a que ha visto el siguiente titular en un famoso diario chileno: “El aumento en la compra de helados tiene una alta correlación con la muerte de personas en Santiago”. ¿Que le recomendaría a su amigue sobre el titular leído?, ¿Debería preocuparse tanto?.
Le recomendaría que debería dejar de preocuparse, dado que el titular divulga probablemente información falsa y sin un estudio estadístico exhastivo.
En la disciplina de la estadística es muy dificil saber cuando tenemos causalidades directas entre variables, dado que varias relaciones ocultan otras realidades y relaciones ocultas, las cuales llamamos variables de confusión.
Lo anterior lleva a que la Correlación entre variables no implica Causalidad, lo cuál en palabras más simples y aplicado al ejemplo significa que no porque aumente la cantidad de ofertas de helados significa que más gente va a morir. No estamos considerando otras variables ocultas, como por ejemplo, probablemente todo lo anterior se debe al aumento de la temperatura en la ciudad de Santiago.
Pregunta 3:
Señale las diferentes aplicaciones que poseen las visualizaciones: Boxplot, histograma, gráfico de pie y scatterplot.
Boxplot: Nos permite visualizar la distribución de datos numéricos, pero a diferencia de los gráficos de funciones de densidad o probabilidad como la normal, este segmenta los datos en cuartiles y es capaz de mostrarlos utilizando una única dimensión. La marca central indica la mediana, el rectángulo central indica donde se ubica el 25% de los datos menores y 25% de los datos mayores a la mediana, mientras que los bigotes indican los valores más externos. Este a su vez sirve para poder observar outliers o valores atípicos dentro de los datos.
Histograma: Nos permiten visualizar la distribución completa de los datos numéricos, ya sean discretos o continuos (utilizando rangos de valores). Al visualizar por completo la distribución de los datos por rangos, podemos tener una idea de la forma que toma la función de de densidad de probabilidad del dataset.
Gráfico de Pie: A diferencia de los gráficos anteriores, el Pie Chart sirve para visualizar proporciones numéricas de un conjunto de datos. Es útil cuando queremos visualizar distribuciones de variables categóricas, pero no es viable para mostrar demasiados datos al dificultar su comprensión.
Scatterplot: Este gráfico utiliza coordenadas cartesianas para graficar datos numéricos con dos o más variables. Es útil para apreciar el nivel de dependencia y correlación entre dos variables numéricas de un dataset, sean lineales o no lineales.
Pregunta 4:
Suponga que esta estudiando la diferencia en los sueldos de las personas que viven en Santiago y Rancagua. Suponiendo que los datos poseen outliers, ¿Que métrica de resumen utilizaría para comparar los datos?. Justifique su respuesta.
Utilizaría la mediana, ya que a diferencia del promedio esta es más robusta ante valores extremos en los datos, debido a que esta depende de la cantidad y densidad de los datos en vez de sus valores.
Pregunta 5:
En base al mismo dataset de sueldos para las regiones de Santiago y Rancagua, le comentan que existe un error en los datos y que estos deben ser modificados aumentando un 10% el valor original y sumando \(15.000\) a cada uno de los datos. ¿Como se ve afectada la media, mediana y desviación estándar con esta modificación?. Explique a través de ecuaciones el cambio que experimentan las métricas de resumen respecto al valor original, considere para el caso de la media \(\bar{X}_{old} = \dfrac{1}{m} \sum^{m}_{i=1} x_i\) y \(sd_{old} = \sqrt{\dfrac{1}{(m-1)}\sum_{i=1}^{m}(x_i-\bar{x})^{2}}\) para la desviación estándar.
Promedio:
\(\bar{X}_{new} = \dfrac{1}{m} \sum^{m}_{i=1} 1.1x_i + 15000 = \dfrac{1.1}{m} \sum^{m}_{i=1} x_i + \dfrac{1}{m} \sum^{m}_{i=1} 15000\)
Donde la sumatoria de los x de al final corresponde justo al promedio anterior, teniendo que:
\(\bar{X}_{new} = 1.1 \bar{X}_{old} + 15000\)
Que es justamente la modificación realizada a cada dato por separado. Luego este cambió linealmente respecto a la modificación.
Mediana: Debido a que esta depende de la cantidad de datos y no de sus valores, esta variará de la misma forma en la que lo hace la mediana anterior. Luego la mediana también cambia linealmente respecto a la modificación.
Desviación Estándar:
\(sd_{new} = \sqrt{\dfrac{1}{(m-1)}\sum_{i=1}^{m}((1.1x_i + 15000) - (1.1\bar{x} + 15000)^{2}} = \sqrt{\dfrac{1.1^{2}}{(m-1)}\sum_{i=1}^{m}(x_i - \bar{x})^{2}}\)
Donde la sumatoria final corresponde justamente a la desviación estandar anterior pero por un factor 1.1, teniendo que:
\(sd_{new} = 1.1 sd_{old}\)
Que es el la pendiente de la modificiación realizada a los datos, Luego este aumentó en un 10 porciento.
Pregunta 6:
Suponga que debe responder un examen sorpresa de 10 preguntas, con 5 alternativas por cada pregunta. ¿Cual es la probabilidad de obtener mas de 5 alternativas correctas si responde de forma aleatoria todo el examen?.
Nota: Puede resolver el ejercicio desarrollándolo a mano o utilizando código en R.
# Debido a que queremos calcular la probabilidad de un suceso dicotómico, es claro
# que tenemos una distribución Binomial:
# La probabilidad de responder cada pregunta correctamente es de un quinto:
p = 0.2
# Sabemos que son más de 5 preguntas de 10:
x = c(6:10)
size = 10
# Luego la probabilidad solicitada es:
prob = 0
for(val in x){
prob = prob + dbinom(val, size, p)
}
prob## [1] 0.006369382Pregunta 7:
Supongamos que el 10% de los alumnos del curso utilizan Macintosh, el 60% utiliza Windows y el 30% utiliza Linux. Supongamos que el 50% de los usuarios de Mac, el 78% de los usuarios de Windows y el 20% de los usuarios de Linux han sucumbido bajo un terrible virus. Al seleccionar una persona al azar nos enteramos de que su sistema está infectado por el virus. ¿Cuál es la probabilidad de que sea un alumno con Windows?.
Recopilemos las probabilidades conocidas:
SO: M 10%, W 60%, 30$ L.
Virus: M 50%, 78% W, 20% L.
Sabiendo que la persona escogida si tiene un virus, debemos calcular la probabilidad de que el alumno tenga windows, esto es, calcular \(\mathbb{P}( \text{Windows} | \text{Virus})\). Ahora por Teorema de Bayes, tenemos que \(\mathbb{P}(X | Y) = \frac{\mathbb{P}(Y | X)*\mathbb{P}(X) }{\mathbb{P}(Y)}\), y definiendo \(X\) como el evento de que alguien tenga Windows, e \(Y\) como el evento de que alguien tenga el Virus, tendremos lo siguiente:
\(\mathbb{P}(\text{Windows} | \text{Virus}) = \frac{\mathbb{P}(\text{Virus} | \text{Windows}) * \mathbb{P}(\text{Windows})}{\mathbb{P}(\text{Virus})}\)
\(\mathbb{P}(\text{Windows} | \text{Virus}) = \frac{0.78 * 0.6}{\mathbb{P}(\text{Virus})}\)
Ahora, usando que:
\(\mathbb{P}(X) = \sum^{n}_{i=1} \mathbb{P}(X, Y_{i}) = \sum^{n}_{i=1} \mathbb{P}(X | Y_{i}) * \mathbb{P}(Y_{i})\)
Podemos calcular \(\mathbb{P}(\text{Virus}) = \sum^{n}_{i=1} \mathbb{P}(X | Y_{i}) * \mathbb{P}(Y_{i}) = 0.5*0.1 + 0.78*0.6 + 0.2*0.3 = 0.578\)
Finalmente \(\mathbb{P}(\text{Windows} | \text{Virus}) = \frac{0.78 * 0.6}{0.578} \approx 0.81\)
Pregunta 8:
Señale si las siguientes declaraciones son verdaderas o falsas respecto a las variables aleatorias:
- Como las variables aleatorias son funciones que nos permiten obtener valores de probabilidad, siempre podemos obtener \(\mathbb{P}(X=x)>0\) evaluando en una \(f(x)\) continua y discreta.
- Una PDF bien definida solo puede tener valores menores a 1 y un área debajo de la curva igual a 1.
- La CDF puede ser representada como la integral de la PDF y PMF.
- Una CDF es definida para todo x, continua hacia la derecha y no es decreciente.
Falso, dado que para una \(f(x)\) continua \(\mathbb{P}(X=x)\) siempre es 0, solo pueden evaluarse en rangos.
Falso, ya que una PDF si puede tener valores iguales o mayores a 1.
Falso, ya que en el caso de las PMF estamos con variables discretas, por lo que no es una función integrable. En vez de eso es necesaria una sumatoria.
Verdadera.
Pregunta 9:
Una famosa fabrica de dulces señala que solo el \(5\%\) de sus dulces contienen menos de \(350\) gramos. Si los dulces elaborados por la fabrica distribuyen de forma normal, con media \(\mu\) y desviación estándar \(11.2\). Responda las siguientes preguntas:
- Encuentre la media del producto.
- Señale el porcentaje de dulces que se encuentran sobre los \(390\) gramos.
Nota: Puede ser útil https://www.statskingdom.com/z_table.html
Utilizando el enlace podemos determinar que \(\mathbb{P}(X < -1.645) \approx 0.05\). Luego utilizando que \(X = \mu + \sigma Z\) tenemos que:
\(350 = \mu - 1.645 \sigma\)
\(\mu = 350 + 1.645*11.2 \approx 368.4\). Finalmente \(\mu \approx 368.4\) (a)
Ahora debemos transformar de la distribución actual a la normal, lo cuál se hace utilizando \(Z = \frac{X - \mu}{\sigma}\). Así desarrollamos como sigue:
\(Z = \frac{390 - 368.4}{11.2} \approx 1.93\). Y viendo la tabla del enlace podemos determinar que \(\mathbb{P}(X < 1.93) \approx 0.97 \Rightarrow \mathbb{P}(X > 1.93) \approx 1 - 0.97 = 3 \%\) (b)
Segunda Parte: Elaboración de Código
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Pregunta 1: Visualización de Datos
Para esta pregunta usted deberá trabajar en base al conjunto de datos hearth_database.csv, el cual esta compuesto por las siguientes variables:
- target: Señala si el paciente tuvo un infarto.
- sex: Sexo de los sujetos de prueba.
- fbs: Azúcar en la sangre con ayunas. Esta variable señala solo si se encuentra <=120 o >120.
- exang: Angina de pecho inducida por el ejercicio.
- cp: Tipo de dolor de pecho.
- restecg: Resultados electrocardiográficos en reposo.
- slope: Pendiente del segmento ST máximo de ejercicio.
- ca: Número de buques principales.
- thal: Thalassemia.
- age: Edad en años.
- trestbps: Presión arterial en reposo.
- chol: colesterol sérico en mg/dl.
- thalach: Frecuencia cardíaca máxima alcanzada.
- oldpeak: Depresión del ST inducida por el ejercicio en relación con el reposo.
En base al dataset propuesto realice un análisis exploratorio de los datos (EDA). Para su análisis enfoquen el desarrollo en las siguientes tareas:
- Obtenga la media, mediana, quintiles y valores máximos desde los datos que componen el dataset.
- Obtenga la Matriz de correlación de Pearson y visualice la relación entre las variables numéricas.
- Visualice los boxplot para las variables numéricas.
- Visualice a través de un histograma como distribuyen las variables respecto a los TARGET.
Respuesta
Primero que todo, leemos el dataset y lo guardamos en una variable para poder trabajar con estos. Para verificar que efectivamente guardamos el dataset podemos utilizar la función head() que entrega la cabeza del dataframe:
hdb <- read.csv("hearth_database.csv")
head(hdb)Es importante chequear si es que existen valores nulos dentro del dataset antes de realizar la exploración sobre este. Para revisar la presencia de valores inexistentes se utilizará la función is.na(), que retorna un valor de verdad dependiendo si el dato es NA o no. Realizaremos la operación suma sobre estos valores booleanos, si la suma resultante es 0 quiere decir que no se encontraron valores en el dataset en que is.na() retornará TRUE, lo que nos indicaría que no hay datos inexistentes en las columnas.
sum(is.na(hdb))## [1] 0Ahora, la función summary aplica estadísticas a cada columna. En particular, indica el promedio, mediana, quantiles, máximo, mínimo, entre otros. Además, cambiaremos el nombre de la columna age que tenía un error de tipeo:
colnames(hdb)[10] <- "age"
summary(hdb)## target sex fbs exang
## Length:303 Length:303 Length:303 Length:303
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## cp restecg slope ca
## Length:303 Length:303 Min. :0.000 Min. :0.0000
## Class :character Class :character 1st Qu.:1.000 1st Qu.:0.0000
## Mode :character Mode :character Median :1.000 Median :0.0000
## Mean :1.399 Mean :0.7294
## 3rd Qu.:2.000 3rd Qu.:1.0000
## Max. :2.000 Max. :4.0000
## thal age trestbps chol
## Min. :0.000 Min. :29.00 Min. : 94.0 Min. :126.0
## 1st Qu.:2.000 1st Qu.:47.50 1st Qu.:120.0 1st Qu.:211.0
## Median :2.000 Median :55.00 Median :130.0 Median :240.0
## Mean :2.314 Mean :54.37 Mean :131.6 Mean :246.3
## 3rd Qu.:3.000 3rd Qu.:61.00 3rd Qu.:140.0 3rd Qu.:274.5
## Max. :3.000 Max. :77.00 Max. :200.0 Max. :564.0
## thalach oldpeak
## Min. : 71.0 Min. :0.00
## 1st Qu.:133.5 1st Qu.:0.00
## Median :153.0 Median :0.80
## Mean :149.6 Mean :1.04
## 3rd Qu.:166.0 3rd Qu.:1.60
## Max. :202.0 Max. :6.20Si de otro modo quisieramos todos los quintiles de los datos numéricos, es necesario aplicar la función quantile() a cada columna del dataframe, lo cual podemos hacer con la conocida función lapply() con las columnas como parámetros y la función quantile como la función a aplicar a cada columna:
quint <- function(x) quantile(x, probs = seq(0, 1, 1/5))
quintiles <- lapply(hdb[7:14], quint)
quintiles## $slope
## 0% 20% 40% 60% 80% 100%
## 0 1 1 2 2 2
##
## $ca
## 0% 20% 40% 60% 80% 100%
## 0 0 0 1 2 4
##
## $thal
## 0% 20% 40% 60% 80% 100%
## 0 2 2 2 3 3
##
## $age
## 0% 20% 40% 60% 80% 100%
## 29 45 53 58 62 77
##
## $trestbps
## 0% 20% 40% 60% 80% 100%
## 94 120 126 134 144 200
##
## $chol
## 0% 20% 40% 60% 80% 100%
## 126.0 204.0 230.0 254.0 285.2 564.0
##
## $thalach
## 0% 20% 40% 60% 80% 100%
## 71 130 146 159 170 202
##
## $oldpeak
## 0% 20% 40% 60% 80% 100%
## 0.00 0.00 0.38 1.12 1.90 6.20Ahora, para visualizar la matriz de correlación de Pearson, basta ocupar la función cor() utilizando el método pearson. Esto lo haremos a los valores numéricos del dataframe y rondearemos los valores obtenidos a solo 3 decimales:
pearson_matrix = cor(hdb[,7:14], method = "pearson")
round(pearson_matrix, 3)## slope ca thal age trestbps chol thalach oldpeak
## slope 1.000 -0.080 -0.105 -0.169 -0.121 -0.004 0.387 -0.578
## ca -0.080 1.000 0.152 0.276 0.101 0.071 -0.213 0.223
## thal -0.105 0.152 1.000 0.068 0.062 0.099 -0.096 0.210
## age -0.169 0.276 0.068 1.000 0.279 0.214 -0.399 0.210
## trestbps -0.121 0.101 0.062 0.279 1.000 0.123 -0.047 0.193
## chol -0.004 0.071 0.099 0.214 0.123 1.000 -0.010 0.054
## thalach 0.387 -0.213 -0.096 -0.399 -0.047 -0.010 1.000 -0.344
## oldpeak -0.578 0.223 0.210 0.210 0.193 0.054 -0.344 1.000Para visualizar la matriz de correlación utlizaremos la librería corrplot que contiene la función corrplot() para visualizar dependencias entre variables numéricas con colores y circulos de distintos tamaños:
library(corrplot)## corrplot 0.90 loadedcorrplot(pearson_matrix, type = "upper", order = "hclust",
tl.col = "black", tl.srt = 50)
Para visualizar los Boxplots de las variables numéricas podemos generar un gráfico simple de boxplot con el rango de las columnas que queremos graficar, obteniendo el siguiente resultado:
boxplot(x=hdb[,7:14],main="Boxplots")
Para poder apreciar de mejor manera las variables con valores bajos, podríamos separar las columnas slope, ca, thal y oldpeak en un único gráfico, debido a que sus rangos son parecidos y no se aprecian en conjunto con el resto. Así tenemos:
boxplot(x=hdb[, c(7,8,9,14)], main="Low Values Boxplots")
boxplot(x=hdb[, c(10,11,12,13)], main="High Values Boxplots")
Ahora, pasaremos a visualizar todas las variables numéricas en función de la varible TARGET a través de histogramas:
library(ggplot2)
library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.0.5library(grid)
plot1 <- ggplot(hdb[c(1,7)], aes(x=slope)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot2 <- ggplot(hdb[c(1,8)], aes(x=ca)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot3 <- ggplot(hdb[c(1,9)], aes(x=thal)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot4 <- ggplot(hdb[c(1,10)], aes(x=age)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot5 <- ggplot(hdb[c(1,11)], aes(x=trestbps)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot6 <- ggplot(hdb[c(1,12)], aes(x=chol)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot7 <- ggplot(hdb[c(1,13)], aes(x=thalach)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
plot8 <- ggplot(hdb[c(1,14)], aes(x=oldpeak)) +
geom_histogram(color="black", fill="white") +
facet_grid(target ~ .)
# Hacemos un grid arrange sobre los 8 gráficos para juntarlos en uno solo
grid.arrange(plot1, plot2, plot3, plot4,
plot5, plot6, plot7, plot8,
# Numero de filas y columnas en los que los queremos sortear
ncol=4, nrow=2,
# Títulos
top = textGrob("Hearth variables sorted by target (Y/N) \n",
gp=gpar(fontsize=15)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Pregunta 2: Teorema Central del Limite
Pruebe el teorema central del limite aplicando un muestreo de la media en las distribuciones Poisson, Exponencial y una a su elección. Grafique los resultados obtenidos y señale aproximadamente el numero de muestreos necesarios para obtener el resultado esperado, pruebe esto con las siguientes cantidades de muestreo \(\{10,100,1000,5000\}\). ¿El efecto ocurre con el mismo número de muestreo para todas las distribuciones?.
Respuesta La idea será realizar iteraciones sobre las medias calculadas por cada nueva distribución:
# Definición de variables o estructuras necesarias para el muestreo.
m <- c(10, 100, 1000, 5000)
main_pois <- c()
main_exp <- c()
main_norm <- c()
for (i in 1:4){
n = m[i]
pois <- c()
exp <- c()
norm <- c()
# Realizar el muestreo de las distribuciones.
for(i in 1:n) {
# Elementos aleatorios con distribución de Poisson
poisson <- rpois(1000, lambda = 3)
# Elementos aleatorios con distribución Exponencial
exponential <- rexp(1000, rate = 1)
# Elementos aleatorios con distribución Normal
normal <- rnorm(1000, 0, 1)
pois <- append(pois, mean(poisson))
exp <- append(exp, mean(exponential))
norm <- append(norm, mean(normal))
}
main_pois <- append(main_pois, pois)
main_exp <- append(main_exp, exp)
main_norm <- append(main_norm, norm)
}Ya teniendo todos nuestros valores aleatorios, pasaremos a graficar las medias en histogramas con el fin de observar la distribución y densidad de la variabilidad de las medias calculadas.
# Primero creamos un dataframe de los datos aleatorios, con atributo 'vals'
pois_dt <- data.frame(vals=main_pois[1:10])
# Posteriormente creamos el plot que mostraremos en pantalla
p1 <- ggplot(pois_dt, aes(x=vals)) + # Elegimos el dataframe a plotear
geom_histogram(aes(y=..density..), colour="black", fill="white") + # Creamos el Histograma
geom_density(alpha=.2, fill="#FF6666") + # Ploteamos la densidad del dataframe
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) + # Ploteamos la mediana del dataframe
geom_vline(aes(xintercept=3),
color="red", linetype="solid", size=1) + # Ploteamos además el lambda
ggtitle("n = 10") +
theme(plot.title = element_text(hjust = 0.5)) # Creamos el título y su tema
pois_dt <- data.frame(vals=main_pois[11:110])
p2 <- ggplot(pois_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=3),
color="red", linetype="solid", size=1) +
ggtitle("n = 100") +
theme(plot.title = element_text(hjust = 0.5))
pois_dt <- data.frame(vals=main_pois[111:1110])
p3 <- ggplot(pois_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=3),
color="red", linetype="solid", size=1) +
ggtitle("n = 1000") +
theme(plot.title = element_text(hjust = 0.5))
pois_dt <- data.frame(vals=main_pois[1111:6110])
p4 <- ggplot(pois_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=3),
color="red", linetype="solid", size=1) +
ggtitle("n = 5000") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2,
top = textGrob("Poisson Distribution \n lambda = 3",
gp=gpar(fontsize=15)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
exp_dt <- data.frame(vals=main_exp[1:10])
e1 <- ggplot(exp_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=1),
color="red", linetype="solid", size=1) +
ggtitle("n = 10") +
theme(plot.title = element_text(hjust = 0.5))
exp_dt <- data.frame(vals=main_exp[11:110])
e2 <- ggplot(exp_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=1),
color="red", linetype="solid", size=1) +
ggtitle("n = 100") +
theme(plot.title = element_text(hjust = 0.5))
exp_dt <- data.frame(vals=main_exp[111:1110])
e3 <- ggplot(exp_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=1),
color="red", linetype="solid", size=1) +
ggtitle("n = 1000") +
theme(plot.title = element_text(hjust = 0.5))
exp_dt <- data.frame(vals=main_exp[1111:6110])
e4 <- ggplot(exp_dt, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=1),
color="red", linetype="solid", size=1) +
ggtitle("n = 5000") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(e1, e2, e3, e4, ncol=2, nrow=2,
top = textGrob("Exponential Distribution \ nrate = 1",
gp=gpar(fontsize=15)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
norm_df <- data.frame(vals=main_norm[1:10])
n1 <- ggplot(norm_df, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=0),
color="red", linetype="solid", size=1) +
ggtitle("n = 10") +
theme(plot.title = element_text(hjust = 0.5))
norm_df <- data.frame(vals=main_norm[11:110])
n2 <- ggplot(norm_df, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=0),
color="red", linetype="solid", size=1) +
ggtitle("n = 100") +
theme(plot.title = element_text(hjust = 0.5))
norm_df <- data.frame(vals=main_norm[111:1110])
n3 <- ggplot(norm_df, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=0),
color="red", linetype="solid", size=1) +
ggtitle("n = 1000") +
theme(plot.title = element_text(hjust = 0.5))
norm_df <- data.frame(vals=main_norm[1111:6110])
n4 <- ggplot(norm_df, aes(x=vals)) +
geom_histogram(aes(y=..density..), colour="black", fill="white") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(vals)),
color="blue", linetype="solid", size=1) +
geom_vline(aes(xintercept=0),
color="red", linetype="solid", size=1) +
ggtitle("n = 5000") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(n1, n2, n3, n4, ncol=2, nrow=2,
top = textGrob("Normal Distribution \n mu = 0, sigma = 1",
gp=gpar(fontsize=15)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Notemos que, pese a que todas presentan la distribución Normal esperada por el Teorema Central del Límite, estas distribuciones tienen formas ligeramente distintas entre sí. Por ejemplo, el gráfico generado por las distribuciónes exponencial y Poisson contienen más valores externos hacia la derecha que hacia la izquierda (con n = 100 y 1000), mientras los generados por la distribución normal no. Esto puede deberse a que la distribución exponencial y la de Poisson tiene valores más externos a la curva principal, por lo que es de esperar que las medias calculadas tiendan a ser levemente mayores al resto mientras la proporción de elementos extremos sea relativamente alta al resto (lo cuál ocurre a menor cantidad de datos). Además, podemos apreciar que mientras más aumenta n, la distribución de medias presenta menor proporción de valores externos y más se parece a la distribución normal teórica, lo cuál corresponde justamente a una implicancia del Teorema Central del Límite.
Pregunta 3: Ley de los Grandes Numeros.
Lanzamiento de monedas
Realice el experimento de lanzar una moneda cargada 1000 veces y observe el comportamiento que tiene la probabilidad de salir cara. Para realizar el experimento considere que la moneda tiene una probabilidad de \(4/5\) de salir cara. Grafique el experimento para las secuencias de intentos que van desde 1 a 1000, señalando el valor en que converge la probabilidad de salir cara.
Respuesta
pcara = c(0)
caras = 0
for (lanzamientos in 1:1000) {
coin <- sample(1:10, 1) # Tomaremos un numero al azar de 1 al 10
if (coin < 9){
caras <- caras + 1
}
prob = caras/lanzamientos
pcara <- append(pcara, prob)
}
probs <- data.frame(probabilities = pcara)
# Gráfico de la convergencia
graph <- ggplot(probs, aes(x=c(0:1000), y=probabilities)) + # graficamos el dataframe
geom_line() + # graficamos los datos en forma de línea
ylim(0, 1) + # ponemos límites a los ejes
xlim(0, 1000) +
geom_hline(aes(yintercept=0.8),
color="red", linetype="solid", size=1, alpha = 0.4) + # Ponemos una línea en 4/5
xlab("Lanzamientos") +
ylab("Probabilidad") +
ggtitle("Probabilidad acumulada de salir cara") +
theme(plot.title = element_text(hjust = 0.5))
graph
El problema de Monty Hall
Remontándonos en la televisión del año 1963, en USA existía un programa de concursos donde los participantes debían escoger entre 3 puertas para ganar un premio soñado. El problema del concurso era que solo detrás de 1 puerta estaba el premio mayor, mientras que detrás de las otras dos habían cabras como “premio”.
Una de las particularidades de este concurso, es que cuando el participante escogía una puerta, el animador del show abría una de las puertas que no fue escogida por el participante (Obviamente la puerta abierta por el animador no contenía el premio). Tras abrir la puerta, el animador consultaba al participante si su elección era definitiva, o si deseaba cambiar la puerta escogida por la otra puerta cerrada.
Imagine que usted es participante del concurso y desea calcular la probabilidad de ganar el gran premio si cambia de puerta en el momento que el animador se lo ofrece. Utilizando listas/arrays/vectores simule las puertas del concurso, dejando aleatoriamente el premio en alguna posición del array. Hecho esto, genere un numero de forma aleatoria para escoger una de las puerta (posiciones de la estructura), para luego ver si cambiando de posición tendrá mayores posibilidades de ganar el premio. Genere N veces el experimento y grafique cada una de las iteraciones, tal como se hizo en el ejercicio de las monedas.
Respuesta:
# Creamos una función que simule el juego
montyhall <- function(cambiar = TRUE){
Puertas <- sample(1:3,3) # Puertas donde la posición que tiene el 3 es el premio
posicion <- sample(1:3,1) # Elección del participante.
if (cambiar){
if (Puertas[posicion] == 3){
Eleccion <- 0
} else {
Eleccion <- 1
}
} else {
if (Puertas[posicion] == 3){
Eleccion <- 1
} else {
Eleccion <- 0
}
}
return(Eleccion) # Retornamos la elección, esta puede que tenga el premio o no
}
# Función que simula N juegos
n_juegos <- function(n = 10 ,cambiar_puerta = TRUE){
probs <- c(0)
wins <- 0
for (i in c(1:n)){
eleccion <- montyhall(cambiar_puerta)
if (eleccion == 1){
wins <- wins + 1
}
prob = wins/i
probs <- append(probs, prob)
}
return(probs) # Vector de largo n + 1
}Una primera impresión es que la probabilidad de éxito después de cambiar la primera puerta escogida sigue siendo la misma que no cambiarla, pero este razonamiento es erroneo. La probabilidad de éxito eligiendo solo una puerta es claramente de 1/3 y esta sigue manteniendose si en la segunda instancia se decide no cambiar la puerta, pero la probabilidad de éxito al cambiar de puerta en la segunda instancia es de 2/3, debido a que la probabilidad de que la primera puerta escogida sea el premio mayor es de 1/3 lo que implica que la probabilidad de que la segunda puerta tenga el premio mayor es de 1 - 1/3 = 2/3.
Con lo anterior, graficaremos la probabilidad acumulada por cada iteración del juego y además la constante 2/3 para observar cómo esta converge a este valor en vez de 1/3.
n <- 1000
Probs <- data.frame(probabilities = n_juegos(n, TRUE)) # tiene n + 1 elementos
# Gráfico de la convergencia
graph <- ggplot(Probs, aes(x=c(0:n), y=probabilities)) + # graficamos el dataframe
geom_line() + # graficamos los datos en forma de línea
ylim(0, 1) + # ponemos límites a los ejes
xlim(0, n) +
geom_hline(aes(yintercept=2/3),
color="red", linetype="solid", size=1, alpha = 0.4) +
geom_hline(aes(yintercept=1/3),
color="blue", linetype="solid", size=1, alpha = 0.4) +
xlab("Número de Juego") +
ylab("Probabilidad") +
ggtitle("Probabilidad acumulada de ganar al cambiar puerta en Monty Hall") +
theme(plot.title = element_text(hjust = 0.5))
graph
Pregunta 4: ¿Independencia?
Ustedes disponen de los dados D1 y D2, los cuales no están cargados y son utilizados para comprobar que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\) cuando el evento A es independiente del B. Para estudiar la independencia considere que los eventos A y B se definen de la siguiente manera; sea A el evento dado por los valores obtenidos en el lanzamiento del dado D1, este está compuesto por \(A=\{D1=1,D1=2,D1=6\}\). Por otro lado, el evento B viene dado por los valores obtenidos con el dado D2, el que está conformado por \(B=\{D2=1,D2=2,D2=3,D2=4\}\). Con esto, tendremos un \(\mathbb{P}(A)=1/2\), \(\mathbb{P}(𝐵)=2/3\) y \(\mathbb{P}(AB)=1/3\). Compruebe de forma gráfica que al realizar 1000 lanzamientos (u otro valor grande que usted desea probar) se visualiza que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\).
Hecho lo anterior, compruebe el comportamiento de un segundo grupo de eventos, dados por el lanzamiento de solo el dado D1. Donde, los eventos para D1 quedan definidos como: \(A =\{D1=1,D1=2,D1=6\}\) y \(B=\{D1=1,D1=2,D1=3\}\). ¿Se observa independencia en este experimento?.
Se le aconseja que para simular los lanzamientos de dados utilice la función sample() para generar valores aleatorios entre 1 y 6. Compruebe los números generados por la función con los casos favorables de cada uno de los eventos a ser estudiados.
Respuesta:
# Primer grupo de eventos
N_lan = 1000 # Numero de lanzamientos
L_A = c(0) # Lanzamientos favorables A = c(1, 2, 6)
L_B = c(0) # Lanzamientos favorables B = c(1, 2, 3, 4)
L_AB = c(0) # Lanzamientos favorables AB = c(1, 2)
# Contadores para los lanzamientos
C_A <- 0
C_B <- 0
C_AB <- 0
# Ciclo for de simulación de los dados.
for (i in c(1:1000)){
D1 <- sample(1:6,1)
D2 <- sample(1:6,1)
if (D1 %in% c(1, 2, 6)){ # Evento A
C_A <- C_A + 1
}
if (D2 %in% c(1, 2, 3, 4)){ # Evento B
C_B <- C_B + 1
}
if (D1 %in% c(1, 2, 6) && D2 %in% c(1, 2, 3, 4)){ # Evento A y B
C_AB <- C_AB + 1
}
# Guardamos en cada iteración las probabilidades
L_A <- append(L_A, C_A/i)
L_B <- append(L_B, C_B/i)
L_AB <- append(L_AB, C_AB/i)
}
# Creamos tres dataframes con distintos ids
probs_A <- data.frame(probabilities = L_A, id = 1)
probs_B <- data.frame(probabilities = L_B, id = 2)
probs_AB <- data.frame(probabilities = L_AB, id = 3)
# Generamos el gráfico con los tres resultados
g <- ggplot(NULL, aes(c(0:1000), y=probabilities)) + # Ploteamos null de 0 a 1000
ylim(0, 1) +
xlim(0, 1000) +
geom_line(data = probs_A, col = "red") + # Ploteamos encima el dataset del evento A
geom_line(data = probs_B, col = "blue") + # Lo mismo con B
geom_line(data = probs_AB, col = "green") + # Lo mismo con A y B
# Además añadimos constantes con los valores de las probabilidades teóricas
geom_hline(aes(yintercept=1/2),
color="red", linetype="solid", size=1, alpha = 0.3) +
geom_hline(aes(yintercept=2/3),
color="blue", linetype="solid", size=1, alpha = 0.3) +
geom_hline(aes(yintercept=1/3),
color="green", linetype="solid", size=1, alpha = 0.3) +
xlab("Número de Lanzamiento") +
ylab("Probabilidad") +
ggtitle("¿Independencia?") +
theme(plot.title = element_text(hjust = 0.5))
g
Claramente podemos apreciar que las probabilidades calculadas calzan con las observadas en el experimento (como consecuencia de la Ley de los Grandes Números), lo cual confirma la independencia de ambos eventos A y B. Por otro lado, utilizando la misma idea de ploteo, veamos el segundo grupo de eventos que solo involucra a un dado.
# Segundo grupo de eventos
N_lan = 1000 # Numero de lanzamientos
L_A = c(0) # Lanzamientos favorables A = c(1, 2, 6)
L_B = c(0) # Lanzamientos favorables B = c(1, 2, 3)
L_AB = c(0) # Lanzamientos favorables AB = c(1, 2)
# Contadores para los lanzamientos
C_A <- 0
C_B <- 0
C_AB <- 0
for (i in c(1:1000)){
D1 <- sample(1:6,1)
if (D1 %in% c(1, 2, 6)){
C_A <- C_A + 1
}
if (D1 %in% c(1, 2, 3)){
C_B <- C_B + 1
}
if (D1 %in% c(1, 2, 6) && D1 %in% c(1, 2, 3)){
C_AB <- C_AB + 1
}
L_A <- append(L_A, C_A/i)
L_B <- append(L_B, C_B/i)
L_AB <- append(L_AB, C_AB/i)
}
probs_A <- data.frame(probabilities = L_A, id = 1)
probs_B <- data.frame(probabilities = L_B, id = 2)
probs_AB <- data.frame(probabilities = L_AB, id = 3)
g <- ggplot(NULL, aes(c(0:1000), y=probabilities)) +
ylim(0, 1) +
xlim(0, 1000) +
geom_line(data = probs_A, col = "red") +
geom_line(data = probs_B, col = "blue") +
geom_line(data = probs_AB, col = "green") +
# Probabilidad teórica del evento A y del evento B por separado (un medio)
geom_hline(aes(yintercept=1/2),
color="red", linetype="solid", size=1, alpha = 0.3) +
# Probabilidad de A y B asumiendo independencia
geom_hline(aes(yintercept=1/4),
color="green", linetype="solid", size=1, alpha = 0.3) +
# Probabilidad real del evento A y B
geom_hline(aes(yintercept=1/3),
color="black", linetype="solid", size=1, alpha = 0.3) +
xlab("Número de Lanzamiento") +
ylab("Probabilidad") +
ggtitle("¿Independencia?") +
theme(plot.title = element_text(hjust = 0.5))
g
Observando las rectas teóricas podemos ver claramente que la probabilidad del evento AB no sigue la probabilidad esperada de 25% (recta verde), si no que ronda es del 33% (según la recta constante de color negro). Esto nos indica que los eventos A y B NO son eventos independientes, lo cuál tiene sentido debido a que ambos eventos dependen del valor tomado por un único dado, donde los valores son claramente excluyentes unos con otros (que salga 2 y 3 por ejemplo son eventos dependientes, pues no es posible que salga 3 si es que salió 2).
Pregunta 5: La Ruina del Jugador
Un amigo ludópata suyo le comenta que el truco de jugar en el casino esta en no parar de apostar y apostando lo mínimo posible. Ya que así, tienes mas probabilidades de ganar el gran pozo que acumula el juego. Usted sabiendo la condición de su amigo, decide no creer en su conjetura y decide probar esto a través de un experimento.
Para realizar el experimento usted decide asumir las siguientes declaraciones, bajo sus observaciones:
- La probabilidad de ganar en un juego del casino es \(9/19\)
- Sabe que su amigo posee fondos en el rango de 0 a 200 dolares.
- Las apuestas como mínimo deben ser igual a 5 dolares.
- El monto de las apuestas no cambia y son siempre igual a la primera. Por ejemplo, si su amigo apuesta 50 dolares, todos los próximos juegos apuesta 50 hasta que se acaba su dinero.
En el experimento deberá obtener la evolución de los fondos hasta que el jugador se queda sin fondos para jugar. Puede ser útil seguir la lógica de una moneda cargada para realizar esto. Pruebe esto con una apuesta igual a 5, 25 y 50 graficando los resultados. Comente los resultados obtenidos.
Respuesta
# Función para obtener el desarrollo de las apuestas
ruina <- function(fondos = 100, apuesta = 5){
vec_fondos <- c(fondos)
while (0<fondos & fondos<200) {
bet <- sample(1:19, 1) # Tomaremos un numero al azar de 1 al 19
if (bet < 10){
fondos <- fondos + apuesta
} else{
fondos <- fondos - apuesta
}
vec_fondos <- append(vec_fondos, fondos)
}
return(vec_fondos) # Devuelve un vector con el desarrollo de los fondos
}
plot(ruina(), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 5)")
plot(ruina(apuesta = 25), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 25)")
plot(ruina(apuesta = 50), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 50)")
Con lo anterior podemos notar que claramente nuestra conjetura es correcta, y que nuestro compañero está equivocado. Realizar siempre la misma apuesta con un valor bajo da la ilusión de ganar más, debido a que esta puede extenderse por más tiempo que con apuestas de mayor valor, además de que se necesita mayor cantidad de aciertos para poder llegar a los límites del fondo.
Por otro lado, con las apuestas de mayor valor se reducen rápidamente la cantidad de apuestas realizadas, debido a que pueden sobrepasar los límites del fondo más rápidamente respecto a las apuestas de menor valor. Por lo anterior podemos inferir que realizar apuestas de valores más altos tienen mayor probabilidad de generar muchas más ganancias por secuencias de victorias consecutivas, pero equivalentemente puede llevar a la quiebra más rápidamente.
A work by CC6104